Configuración Inicial
%%html
<style>
.qst {background-color: #b1cee3; padding:10px; border-radius: 5px; border: solid 2px #5D8AA8;}
.qst:before {font-weight: bold; content:"Pregunta"; display: block; margin: 0px 10px 10px 10px;}
h1, h2, h3 {color: #5D8AA8;}
.text_cell_render p {text-align: justify; text-justify: inter-word;}
</style>
%matplotlib inline
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import colors
from mpl_toolkits.mplot3d import Axes3D
# My imports
import utils
from typing import Callable, Optional, Union, Type
from sklearn.utils.extmath import svd_flip
from sklearn.gaussian_process.kernels import RBF as rbf_kernel
from sklearn.metrics import pairwise_distances
from sklearn.metrics.pairwise import laplacian_kernel
from sklearn.datasets import make_s_curve, make_swiss_roll
from sklearn.neighbors import NearestNeighbors
from sklearn.cluster import KMeans
from sklearn.cluster import SpectralClustering as SpectralClusteringSklearn
from scipy.sparse.csgraph import laplacian as scipy_laplacian
# Import and initialize scaler
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
matplotlib.rc('figure', figsize=(15, 5))
my_cmap = plt.cm.Spectral
seed = 123
my_colors = np.array(['tomato', 'skyblue', 'gold', 'lime', 'orchid'])
Este notebook contiene el enunciado de la segunda parte del examen de Métodos Funcionales en Aprendizaje Automático.
Por favor, lee con cuidado cada uno de los tres ejercicios así como las preguntas que en ellos se formulan, y explica de manera concisa cada asunción y conclusión que hagas.
Debéis entregar un informe (que puede realizarse directamente sobre el notebook) y el código implementado para completar los ejercicios, así como las referencias consultadas en caso de necesidad. También se debe entregar el código de honor firmado. Todo esto se enviará vía Moodle en un .zip.
Solución.-
Asumimos que el peso de los lazos (lados conectado un nodo consigo mismo) es 0, al no aparecer ningún peso en la imagen anterior. Definimos la matriz $W$ y computamos el laplaciano del grafo:
def evaluate_laplacian(W, f):
# Compute diagonal matrix
D = np.diag(np.sum(W, axis=1))
# Compute laplacian
L = D - W
# Evaluate and also return the laplacian
return L, f @ L @ f.T
W = np.array([
[0, 20, 1, 0, 0],
[20, 0, 10, 0, 0],
[1, 10, 0, 0, 0],
[0, 0, 0, 0, 5],
[0, 0, 0, 5, 0]
], dtype=float)
f = np.array([1, 2, 6, 5, 4], dtype=float)
L, result = evaluate_laplacian(W, f)
print('L: {} \n'.format(L))
print('f^T L f: {}'.format(result))
L: [[ 21. -20. -1. 0. 0.] [-20. 30. -10. 0. 0.] [ -1. -10. 11. 0. 0.] [ 0. 0. 0. 5. -5.] [ 0. 0. 0. -5. 5.]] f^T L f: 210.0
Solución.-
W = np.array([
[0, 1, 20, 0, 0],
[1, 0, 10, 0, 0],
[20, 10, 0, 0, 0],
[0, 0, 0, 0, 5],
[0, 0, 0, 5, 0]
], dtype=float)
L, result = evaluate_laplacian(W, f)
print('L: {} \n'.format(L))
print('f^T L f: {}'.format(result))
L: [[ 21. -1. -20. 0. 0.] [ -1. 11. -10. 0. 0.] [-20. -10. 30. 0. 0.] [ 0. 0. 0. 5. -5.] [ 0. 0. 0. -5. 5.]] f^T L f: 666.0
Como podemos ver, el resultado ha cambiado de 210 a 666. Para responder a por qué ocurre este hecho hemos de dar contexto a las operaciones que estamos realizando. En el contexto de los Laplacian Eigenmaps, se crea un grafo para un conjunto de datos con pesos asociados a la afinidad de los nodos. Es decir, cuanto más cerca estén dos muestras $x_i, x_j$, mayor será el peso asociado entre las mismas, $w_{ij}$.
En las notas de clase definimos el error de reconstrucción:
$$ \mathcal J(y) = \frac{1}{2} \sum_{i,j}(y_i - y_j)^2 w_{ij}. $$Para minimizar este error hemos de asignar valores de $y_i$ parecidos a elementos con altos pesos (es decir, si los puntos son cercanos, asignamos valores de la función parecidos para que la reconstrucción tenga sentido). También demostramos en clase que el error de reconstrucción coincide con la evaluación de la siguiente forma cuadrática:
$$ \mathcal J(y) = y^T L y. $$En nuestro caso particular, los valores de $f(x_1)$ y $f(x_2)$ son parecidos ($1$ y $2$ respectivamente), mientras que el valor de $f(x_3) = 6$ dista más de ellos. Por lo tanto, para minimizar el error de reconstrucción, los pesos han de cumplir $w_{1,2} \gg w_{2,3} > w_{1,3}$.
Inicialmente: $w_{1,2} = 20, w_{1,3} = 1$, lo que tiene sentido (los nodos 1 y 2 estarán muy cerca mientras que los nodos 1 y 3 estarán muy lejos). En el segundo caso $w_{1,2} = 1, w_{1,3} = 20$, representando a los nodos 1 y 2 muy lejanos y los nodos 1 y 3 muy cercanos, lo que no encaja para nada con los valores de $f(x)$. Es por ello que el valor de $f^T L f$ (el error de reconstrucción) aumenta considerablemente.
Solución.-
Nota inicial: Realizaremos la minimización para ambas matrices $W$ definidas en los apartados anteriores. La solución obtenida es válida considerando cualquiera de las dos.
Sabemos que el laplaciano del grafo no normalizado ($L = D - W$) es una matriz semidefinida positiva. Esto es, $x^T L x \ge 0$ para todo vector $x$. El vector $f_1 = 0$ obtiene como resultado $f_1^T L f_1 = 0$, minimizando la función.
f_1 = np.zeros(5)
_, result = evaluate_laplacian(W, f_1)
print('f^T L f: {}'.format(result))
f^T L f: 0.0
Puesto que nuestra función $g(x) = x^T L x$ es una forma cuadrática, sabemos que es analítica (en particular dos veces derivable), y que una función dos veces derivable es convexa si y sólo si su matriz Hessiana ($L$, por ser una forma cuadrática) es semidefinida positiva, y es estrictamente convexa si y sólo si $L$ es definida positiva.
Si $f$ fuese estrictamente convexa, tendríamos un único mínimo. Sin embargo, sabemos que $\lambda_0 = 0$ siempre es valor propio de $L$, como vimos en clase. En nuestro caso particular:
np.linalg.svd(L)[1]
array([4.74620776e+01, 1.45379224e+01, 1.00000000e+01, 1.17995953e-15,
0.00000000e+00])
Al tener un valor propio $0$ sabemos que $L$ es semidefinida positiva, luego $g$ no es estrictamente convexa y puede tener otro valor mínimo. Para buscarlo intentamos en primer lugar vectores con valores no nulos en una única posición, pero esto nos devuelve valores de la diagonal multiplicados por un número positivo, por lo que nunca obtendremos un cero.
Hemos de hacer uso de los valores negativos, no situados en la diagonal. Mirando la submatriz inferior derecha de $L$ vemos que si utilizamos el vector $f_2 = (0, 0, 0, 1, 1)$ obtendremos un valor nulo, ya las dos últimas filas (y por tanto las dos últimas columnas) suman $0$ cada una:
f_2 = np.array([0, 0, 0, 1, 1])
_, result = evaluate_laplacian(W, f_2)
print('f^T L f: {}'.format(result))
f^T L f: 0.0
Sin embargo, con esto no basta para asegurar que existan soluciones linealmente independientes, pues $f_1$ y $f_2$ son linealmente dependientes: $f_1 = 0 \cdot f_2$. Mirando la submatriz superior izquierda del laplaciano sin normalizar con la idea anterior en mente vemos que la suma de las tres primeras filas (así como la de las tres primeras columnas) es nula, por lo que también nos sirve la solución $f_3 = (1, 1, 1, 0, 0)$:
f_3 = np.array([1, 1, 1, 0, 0])
_, result = evaluate_laplacian(W, f_3)
print('f^T L f: {}'.format(result))
f^T L f: 0.0
Con esto hemos encontrado dos mínimos de $g$ linealmente independientes: $f_2$ y $f_3$.
Un algoritmo muy utilizado tanto en el ámbito de manifold learning, como en clustering o en análisis de grafos (por ejemplo para detectar comunidades) es el método de Spectral Clustering. Durante las clases del curso no nos ha dado tiempo a profundizar en él, por lo que vamos a utilizarlo para resolver un problema.
Solucion.-
El algoritmo de Spectral Clustering está compuesto por dos pasos clave:
X aplicandoles Laplacian Eigenmaps (LE) para reducir a un espacio de n_components dimensiones. La siguiente celda implementa dicho algoritmo de la forma más genérica posible y utilizando como base la clase de DiffusionMap creada para la práctica. Esta incluye:
La matriz de pesos del grafo asociado, W, (también denominada matriz de afinidad) puede calcularse de distintas formas, elegida por el parámetro affinity:
(Por defecto) Si affinity == 'kernel', se utilizará un kernel para este computo. La posición $w_{ij}$ de la matriz será kernel(x_i, x_j). El kernel utilizado se define a partir del parámetro kernel, pudiendo ser rbf (por defecto), laplacian o una función kernel pasada por parámetro. La longitud del kernel ($\sigma$ o $\gamma$ dependiendo del kernel utilizado) se define a partir del parámetro sigma. Este sigma puede calcularse automáticamente utilizando los valores de median o maximum, como se explico en la práctica de DiffussionMap.
Si affinity == knn, se utilizará el algoritmo de $k$ vecinos más cercanos (KNN) para computar la matriz de afinidad. En particular $w_{ij}$ será $1$ si el elemento $x_j$ está entre los $k$ vecinos más cercano del elemento $x_i$. Como está relación no tiene por qué ser simétrica, tras el cálculo de la $W$ inicial la simetrizamos haciendo la media con su traspuesta. El número de vecinos queda definido por el parámetro n_neighbors.
Si affinity == precomputed, la matriz de afinidad se ha calculado previamente y se suministra a través del parámetro precomputed_W.
Una vez la matriz de pesos ha sido calculada, el método de obtención del Laplaciano (fijado con el parámetro laplacian_alg) varía entre:
laplacian_alg == unnormalized: utilizamos el laplaciano sin normalizar: $L = D - W$.laplacian_alg == random walk: utilizamos el laplaciano normalizado por la izquierda, también denominado random walk: $L = I - D^{-1} W$.laplacian_alg == normalized: utilizamos el laplaciano normalizado: $L = I - D^{-\frac{1}{2}} W D^{-\frac{1}{2}}$.Una vez hemos computado el laplaciano, lo diagonalizamos y proyectamos a un espacio con n_components dimensiones. Para ello nos quedaremos con los n_components vectores propios de $L$ con menor valor propio asociado. Si drop_first == True, no consideraremos el primero vector propio, por estar asociado al valor propio $0$.
Una vez obtenido el embedding, aplicamos el algoritmo K-means de Sklearn de forma directa, donde el número de clusters queda definido por el parámetro n_clusters.
Como nota final, todos los cálculos potencialmente aleatorios tienen la semilla fijada. Esta puede controlarse con el parámetro random_state.
Finalmente, para diseñar e implementar esta clase se han consultado las siguientes referencias:
[1]: Wikipedia - Spectral clustering.
[2]: Sklearn - Spectral clustering.
from sklearn.base import BaseEstimator, TransformerMixin
class SpectralClustering(TransformerMixin, BaseEstimator):
"""
Spectral clustering.
"""
def __init__(self,
n_components: int=8,
n_clusters: int=8,
affinity: str='auto',
n_neighbors: int=10,
kernel: Union[str, Callable[[
np.ndarray, np.ndarray], np.ndarray]] = 'auto',
sigma: Union[str, float] = 'median',
percentile: float=0.5,
drop_first: bool=True,
random_state: int=123,
laplacian_alg: str='random walk',
precomputed_W: np.ndarray=None):
"""
Initializes the Spectral clustering object
----------
kernel : kernel function to compute affinity matrix.
Can be 'rbf', 'laplacian' or a callable
sigma : Sigma to apply to the kernel
n_components: number of components to keep after projection
percentile: percentile (0 to 100) used to compute sigma if needed
Returns
-------
self : returns an instance of self.
"""
# Assignment of the hyper-parameters
self.n_components = n_components
self.n_clusters = n_clusters
self.affinity = affinity
self._kernel = kernel
self.sigma = sigma
self.percentile = percentile
self.n_neighbors = n_neighbors
self._drop_first = drop_first
self.random_state = random_state
self.laplacian_alg = laplacian_alg
self.precomputed_W = precomputed_W
# Raise variables
self._fitted = False
self._fixed_sigma = False
def _determine_sigma(self):
"""
Converts string sigma values to actual float sigma values
"""
if self.sigma == 'median':
self.sigma = np.percentile(pairwise_distances(self.X), 50)
elif self.sigma == 'maximum':
self.sigma = np.percentile(pairwise_distances(self.X), 100)
elif self.sigma == 'percentile':
self.sigma = np.percentile(
pairwise_distances(self.X), self.percentile)
self._fixed_sigma = True
def _determine_kernel(self):
"""
Converts kernel value to final kernel
"""
if not self._fixed_sigma:
raise ValueError("Sigma has not been fixed yet")
if self._kernel in ['auto', 'rbf']:
# Recall that this rbf_kernel is the one from
# Gaussian procesess, which receives sigma as parameter
self._kernel = rbf_kernel(self.sigma)
elif self._kernel == 'laplacian':
self._kernel = lambda X, y = None: laplacian_kernel(
X, y, gamma=1/(2*self.sigma**2))
def _compute_weights_matrix(self):
if self.affinity in ['auto', 'kernel']:
# Determine sigma and kernel
self._determine_sigma()
self._determine_kernel()
# Instantiate the kernel and compute the affinity matrix.
# This will be our weights matrix
return self._kernel(self.X, self.X)
elif self.affinity == 'knn':
nn = NearestNeighbors(n_neighbors=self.n_neighbors).fit(self.X)
W = nn.kneighbors_graph(self.X).toarray()
return (W + W.T) / 2.0
elif self.affinity == 'precomputed':
if self.precomputed_W is None:
raise ValueError('precomputed W can\'t be None if affinity is precomputed.')
return self.precomputed_W
else:
raise ValueError('affinity parameter must be one of the following: ' + \
'\t [auto, kernel, knn]')
def _compute_laplacian(self, W):
# Compute diagonal matrix
D = np.sum(W, axis=1)
n = D.shape[0]
# Invert the diagonal. Since it's a vector containing only the
# diagonal, it's O(n)
D_inv = 1 / D
# Compute the Laplacian following the random walk
# or left normalization method.
if self.laplacian_alg == 'unnormalized':
return D - W
elif self.laplacian_alg == 'random walk':
return np.identity(n) - np.diag(D_inv) @ W
elif self.laplacian_alg == 'normalized':
return np.identity(n) - np.diag(np.sqrt(D_inv)) @ W @ np.diag(np.sqrt(D_inv))
def _compute_embedding(self, X=None):
"""
This function computes the embedding using LE.
"""
if X is None and self.affinity != 'precomputed':
print('X can\'t be None if W is not precomputed.')
# Save data
self.X = X
# Compute the weights matrix
W = self._compute_weights_matrix()
# Compute the laplacian
self.L = self._compute_laplacian(W)
# Since the Laplacian is already symmetric (D and W are so),
# we may directly compute its SVD decomposition.
# Remember that eigenvecs has the eigen-vectors as its columns!
eigenvecs, eigenvals, V = np.linalg.svd(self.L, full_matrices=True)
# Make the eigenvector signs match the ones obtained from Sklearn
eigenvecs, _ = svd_flip(eigenvecs, V)
# Invert the order of the eigenvectors (that is, traverse the columns in reverse order).
# We are interested in the ones associated with lower associated eigenvalues.
eigenvecs = np.flip(eigenvecs, axis=1)
eigenvals = np.flip(eigenvals, axis=0)
# Remove first component, with eigenvalue 1
# Recall that we are assuming that the eigenvals/vecs are sorted,
# they are given as such by the np.linalg.svd
# n_components + 1 since we are already not considering 1st position
if self._drop_first:
embedding = eigenvecs[:, 1:self.n_components+1]
else:
embedding = eigenvecs[:, :self.n_components]
return embedding
def fit(self, X=None, y=None):
"""
Compute the embedding vectors for data X
Parameters
----------
X : array-like of shape [n_samples, n_features]
training set.
y : Ignored
Returns
-------
self : returns an instance of self.
"""
# Compute the embedding usnig LE
self.embedding = self._compute_embedding(X)
# Instantiate and fit K-means object
self.k_means = KMeans(n_clusters=self.n_clusters,
random_state=self.random_state)
self.k_means.fit(self.embedding)
# Update fitted state
self._fitted = True
return self
def fit_transform(self, X=None, y=None):
"""Compute the embedding vectors for data X and transform X.
Parameters
----------
X : array-like of shape [n_samples, n_features]
training set.
y : Ignored
Returns
-------
X_red : array-like, shape (n_samples, n_components)
"""
self.fit(X)
return self.k_means.predict(self.embedding)
En primer lugar, Spectral Clustering es un algoritmo de clustering, no de embedding. Podemos comparar Diffusion Maps con el algoritmo que utiliza Spectral Clustering para obtener el embedding. Esto es, Laplacian Eigenmaps con distintos métodos para crear el laplaciano.
Por un lado, ambos métodos buscan utilizar la geometría intríseca de los datos para reducir la dimensionalidad de los mismos, de forma no supervisada. Sin embargo, la motivacición de los métodos es originalmente completamente distinta: Diffusion maps busca simular un proceso de expansión del calor para estudiar dicha geometría, mientras que Laplacian Eigemaps construye un grafo y busca minimizar el error de reconstrucción, por lo que transforma el problema en minimizar un Lagraniano, aunque llegue a una solución similar en última instancia.
Respecto a cómo llegar a dicha solución, los algoritmos difieren entre si en los siguientes detalles:
La forma de computar la matriz de afinidad es distinta en cada caso. En Difussion Maps utilizábamos $D^{-1}W$ mientras que en nuestro caso (al menos siguiente el método pedido) utilizamos $I - D^{-1}W$.
En Diffusion Maps podemos añadir normalización de densidad utilizando el parámetro $\alpha$ implementado en la práctica. Sin embargo, en Laplacian Eigenmaps esto no es una posibilidad.
Tras obtener los vectores propios de la matriz de afinidad, $(\phi_1, \ldots, \phi_n)$, los algoritmos también se diferencian en cómo obtener las nuevas coordenadas de los puntos. En Laplacian Eigenmaps, el punto $i-ésimo$ tendrá por coordenadas los elementos $i$-ésimos de los vectores propios: $x_i = (\phi_1[i], \ldots, \phi_n[i]) \in \mathbb R^n$. En Difussion maps, estas coordenadas se multiplican por el respectivo valor propio elevado al número de pasos simulados $T$: $x_i = (\lambda_1^T \phi_1[i], \ldots, \lambda_n^T \phi_n[i]) \in \mathbb R^n$.
Finalmente, los vectores propios se ordenan dependiendo de los valores propios: en Laplacian Eigemaps en sentido ascendente mientras que en Diffusion Maps, en sentido descendente. Además, en ambos casos se elimina el primer vector propio antes de devolver las coordenadas (correspondientes a los valores propios $0$ y $1$ respectivamente).
Respecto a los parámetros de cada algoritmo, ambos necesitan la dimensionalidad a la que proyectar, así como los detalles de cómo computar la matriz de afinidad (el kernel utilizado si es que se utiliza un kernel, así como los parámetros del mismo). Sin embargo, Difussion Maps tiene un parámetro adicional, el número de pasos a simular. Esto lo hace algo más complicado de ajustar.
Para responder a las preguntas vamos a hacer algunos experimentos con el dataset proporcionado.
Cargamos el dataset proporcionado
data = np.loadtxt('input/data.txt')
X = data[:,:-1]
y = data[:,-1]
# Preparamos las etiquetas
y[y != 1] = -1
En primer lugar mostramos únicamente el embedding obtenido utilizando Laplacian Eigenmaps proyectado a dos dimensiones, asi como las dos primeras dimensiones (de doce) de los datos originales. Utilizamos el kernel rbf con sigma la mediana para obtener el embedding.
sc = SpectralClustering(
n_components = 2,
affinity = 'kernel',
kernel = 'rbf',
sigma = 'median'
)
X_transformed = sc._compute_embedding(X)
utils.scatter_2D(X, y, title='Original data')
utils.scatter_2D(X_transformed, y, title='Transformed data')
Como podemos apreciar, los datos se separan en algo parecido a dos medias lunas intersecadas. Como hemos estudiado amplicamente, este dataset es relativamente sencillo de separar utilizando algún método no lineal no tirivial (como una SVM con kernel no lineal). Sin embargo, es obvio que a la hora de hacer clustering no obtendremos clusters que encaja particularmente bien con las clases originales.
Volvemos a ejecutar este mismo ejemplo utilizando el algoritmo de KNN para computar la matriz de afinidad con 10 vecinos (el valor por defecto).
sc = SpectralClustering(
n_components = 2,
affinity = 'knn',
)
X_transformed = sc._compute_embedding(X)
utils.scatter_2D(X, y, title='Original data')
utils.scatter_2D(X_transformed, y, title='Transformed data')
En esta ocasión los puntos obtenidos son algo más separables linealmente, aunque la intersección no es trivial y hay varios puntos azules dentro del clúster rojo. Intuitivamente este método funcionará mejor a la hora de hacer clustering.
Motivados por la segunda pregunta del enunciado, realizamos un 'mini-gridsearch' sobre valores de $\sigma$ y $k$ para los métodos anteriores, elegidos a mano. Mostramos a continuación los resultados:
sigmas = [0.5, 0.7, 0.9, 1, 1.2, 1.4, 1.5, 2, 5, 10, 'median', 'maximum']
fig, axes = plt.subplots(3, 4, figsize=(20, 15))
for i, sigma in enumerate(sigmas):
sc = SpectralClustering(
n_components = 2,
kernel = 'rbf',
sigma = sigma,
)
X_transformed = sc._compute_embedding(X)
axes[i//4][i % 4].scatter(X_transformed[:, 0], X_transformed[:, 1],
c=y, cmap=my_cmap)
axes[i//4][i % 4].set_title("RBF SC - $\sigma = ${}".format(str(sc.sigma)))
axes[i//4][i % 4].set_yticklabels([])
axes[i//4][i % 4].set_xticklabels([])
ks = [2, 3, 5, 10, 15, 20, 50, 100]
fig, axes = plt.subplots(2, 4, figsize=(20, 15))
for i, k in enumerate(ks):
sc = SpectralClustering(
n_components = 2,
affinity = 'knn',
n_neighbors = k,
)
X_transformed = sc._compute_embedding(X)
axes[i//4][i % 4].scatter(X_transformed[:, 0], X_transformed[:, 1],
c=y, cmap=my_cmap)
axes[i//4][i % 4].set_title("KNN SC - $k = ${}".format(str(sc.n_neighbors)))
axes[i//4][i % 4].set_yticklabels([])
axes[i//4][i % 4].set_xticklabels([])
Para el primer experimento sobre posibles valores de $\sigma$, vemos como a partir de cierto valor de sigma ($\sim 1$) el embedding se estabiliza y ya no varía. Las dos útlimas gráficas son las del la mediana y el máximo respectivamente. Como se puede apreciar, el embedding obtenido por los valores por defecto es esencialmente el mejor.
Respecto al segundo experimento, vemos como el embedding no llega a estabilizarse para los valores de $k$ probados. Sin embargo, conforme se aumenta lo suficiente, los clusters empiezan a mezclarse significativamente. Se puede apreciar en los valores de $k$ mayores de $10$. Sin embargo, el caso de $k=50$ separa relativamente los clusters, aunque tampoco significativamente mejor que en el caso de RBF.
Como sabemos, los algoritmos de manifold learning tienen como objetivo aprovechar la geometría original de las variedades diferenciales en las que se encuentran los datos para desenrollarlas en menor dimensión. Por lo que estamos observando, no parece que los datos estén contenidos en una variedad diferencial desenrollable en dos dimensiones. Probemos con tres:
sc = SpectralClustering(
n_components = 3,
affinity = 'kernel',
)
X_transformed = sc._compute_embedding(X)
utils.scatter_3D(X, y, title="Original swiss-roll data.")
Tampoco parecen proyectables en una variedad diferenciable tres-dimensional.
A continuación probamos el algoritmo completo de Spectral Clustering, que incluye el paso de obtención de clusters tras computar el embedding.
sc = SpectralClustering(
n_components = 2,
n_clusters = 2,
affinity = 'kernel',
)
y_pred = sc.fit_transform(X)
utils.plot_decision_boundary(sc.embedding, y, sc.k_means, title='Clusters with original labels')
utils.plot_decision_boundary(sc.embedding, y_pred, sc.k_means, title='Clusters with predicted labels (RBF)')
sc = SpectralClustering(
n_components = 2,
n_clusters = 2,
affinity = 'knn',
)
y_pred = sc.fit_transform(X)
utils.plot_decision_boundary(sc.embedding, y, sc.k_means, title='Clusters with original labels')
utils.plot_decision_boundary(sc.embedding, y_pred, sc.k_means, title='Clusters with predicted labels (KNN)')
Como ya predijimos, el clustering no separa las clases en clusters, ya que estaos no separan particularmente bien los datos. Si quisiesemos clasificar los datos deberíamos aplicar alguna técnica de aprendizaje supervisado no lineal, que debería poder separar las clases de forma mucho más sencilla.
Respondemos a las preguntas enunciadas al principio de este ejercicio.
El embedding nos ayuda a separar las clases mejor que si no lo utilizásemos, pero desde luego no separa los datos de forma obvia.
En primer lugar he probado con los parámetros por defecto, tanto los de Sklearn ($\sigma = 1 / X.shape[1]$) como los computados en la práctica anterior (el máximo y la mediana). De la misma forma, para el método de KNN comenzamos utilizando el valor por defecto de Sklearn, $k=10$.
Viendo que estos parámetros no separan claramente las clases, he realizado un grid-search manual para estudiar otros posibles valores de $\sigma$ y $k$, sin encontrar mejores resultados significativos.
# Create S-curve data
N = 1000
X, color = make_s_curve(n_samples=N, random_state=seed)
utils.scatter_3D(X, color, title="Original swiss-roll data.")
A continuación utilizaremos el algoritmo de Spectral Clustering implementado por Sklearn y el nuestro sobre los mismos datos para ver si obtenemos el mismo resultado, tanto para el embedding obtenido como los clusters.
sc_sklearn = SpectralClusteringSklearn(
n_components=2,
random_state=seed,
gamma=1
)
y_pred_sklearn = sc_sklearn.fit_predict(X)
print('Sklearn\' s gamma: ', sc_sklearn.gamma)
sigma = 1.0 / np.sqrt(2 * sc_sklearn.gamma)
print('Sklearn\' s sigma: ', 1.0 / np.sqrt(2 * sc_sklearn.gamma))
# Manually compute the embedding
L = sc_sklearn.affinity_matrix_
lamb, v = np.linalg.eig(L)
v, _ = svd_flip(v, v.T)
id_lamb = lamb.argsort()[::-1]
v_ord = v[:, id_lamb]
X_skl = v_ord[:, 1:3]
utils.scatter_2D(X_skl, color, equal_axis=True,
title='Sklearn\'s embedding - $\sigma = ${}'.format(sigma))
# Our implementation
sc = SpectralClustering(
n_components = 2,
affinity = 'kernel',
kernel = 'rbf',
sigma = sigma,
random_state=seed,
)
X_transformed = sc._compute_embedding(X)
y_pred = sc.fit_transform(X)
utils.scatter_2D(X_transformed, color, equal_axis=True,
title='Our embedding - $\sigma = ${}'.format(sigma))
# Plot 3d S-curves with new clusters
utils.scatter_3D(X, y_pred_sklearn, title='Sklearn\'s clusters - $\sigma = ${}'.format(sigma))
utils.scatter_3D(X, y_pred, title='Our clusters - $\sigma = ${}'.format(sigma))
Sklearn' s gamma: 1 Sklearn' s sigma: 0.7071067811865475
c:\users\jose\appdata\local\programs\python\python37\lib\site-packages\matplotlib\collections.py:206: ComplexWarning: Casting complex values to real discards the imaginary part offsets = np.asanyarray(offsets, float)
En primer lugar cabe destacar que la expresión del kernel RBF y la nuestra difieren en la colocación del parámetro:
$$ f_{sklearn}(x,y) = e^{- \gamma \| x - y\|^2} \quad \quad f_{ours}(x,y) = e^{- \frac{1}{2\sigma^2} \| x - y\|^2} $$Utilizando el valor de parámetro por defecto de Sklearn ($\gamma = 1$) hemos de utilizar $\sigma = \frac{1}{\sqrt 2}$.
Vemos como los embeddings obtenidos son muy parecidos pero no son ambos buenos y realmente parecidos, pero no encajan exactamente. Comentaremos este hecho más adelante.
Respecto a los clusters obtenidos, ambos algoritmos obtienen resultados realmente buenos: los clusters corresponden a secciones de la superficie desenrollada, y no a clusters obtenidos utilizando una distancia euclídea en el espacio original tridimensional. Este era el resultado esperado por nuestro algoritmo y confirma que la implementación es correcta.
En el caso del clustering realizado por Sklearn vemos como en los límites de la S, los puntos se asignan a un cluster distinto, no al más cercano en la superficie. Esto cobra sentido mirando el embedding obtenido, donde los límites de la S quedan relativamente cercanos a la zona central. Si colocamos un centro de cluster en esta zona central, los puntos del límite se asignarán a este clúster, como ocurre en nuestro caso.
Finalmente, comentemos por qué no obtenemos exactamente el mismo embedding (y por tanto el mismo clustering) respecto al obtenido por Sklearn. Tras mucho estudio sobre este hecho, la conclusión final es que Sklearn no implementa exactamente el mismo algoritmo de clustering que nosotros. Para aclarar este hecho he tenido que bucear en la implementación de Sklearn:
spectral_clustering utiliza spectral_embedding para obtener el embedding, como podemos ver en esta línea de la implementación.
spectral_embedding obtiene el laplaciano utilizando la función scipy.spare.csgraph.laplacian de Scipy como podemos ver en esta línea de la implementación. He realizado experimentos en los que sustituía nuestro computo de la matriz $L$ por una llamada a esta función, sin variar nuestro resultado. Por lo tanto, nuestra diferencia no está en el cómputo de dicha matriz.
Tras obtener el laplaciano, en sklearn lo modifican, obtieniendo algo distinto a lo que nosotros buscamos, como podemos ver en esta línea de la implementación. Es por esto que no obtenemos el mismo resultado, ellos no utilizan el mismo laplaciano.
Adicionalmente, cabe destacar que para obtener el embedding realizado por el algoritmo de spectral clustering de Sklearn hemos de obtener a mano su matriz de afinidad y los vectores propios. Sin embargo, no hemos de utilizar los asociados a los menores valores propios sino aquellos asociados a los mayores valores propios, lo que ratifica que la implementación del embedding es distinta.
A continuación realizamos un experimento similar al realizado para el dataset proporcionado, mostrando distintos embeddings para distintos valores de los parámetros, tanto para $\sigma$ como para $k$, para observar si podemos mejorar el embedding obtenido.
sigmas = [0.1, 0.2, 0.3, 0.5, 0.6, 1.0/np.sqrt(2),
0.8, 0.9, 1.0, 5.0, 'median', 'maximum']
fig, axes = plt.subplots(3, 4, figsize=(20, 15))
for i, sigma in enumerate(sigmas):
sc = SpectralClustering(
n_components = 2,
kernel = 'rbf',
sigma = sigma,
)
X_transformed = sc._compute_embedding(X)
axes[i//4][i % 4].scatter(X_transformed[:, 0], X_transformed[:, 1],
c=color, cmap=my_cmap)
axes[i//4][i % 4].set_title("RBF SC - $\sigma = ${}".format(str(sc.sigma)))
axes[i//4][i % 4].set_yticklabels([])
axes[i//4][i % 4].set_xticklabels([])
n_neighbors = [10, 20, 50, 100, 150, 200, 500, 600]
fig, axes = plt.subplots(2, 4, figsize=(20, 15))
for i, nn in enumerate(n_neighbors):
sc = SpectralClustering(
n_components = 2,
affinity = 'knn',
n_neighbors = nn,
)
X_transformed = sc._compute_embedding(X)
axes[i//4][i % 4].scatter(X_transformed[:, 0], X_transformed[:, 1],
c=color, cmap=my_cmap)
axes[i//4][i % 4].set_title("KNN SC - $k = ${}".format(str(sc.n_neighbors)))
axes[i//4][i % 4].set_yticklabels([])
axes[i//4][i % 4].set_xticklabels([])
Por un lado, utilizando el kernel RBF para computar la matriz de afinidad obtenemos formas de $V$ en nuestro embedding que despliegan relativamente bien nuestra superficie (sin desplegarlo en forma de alfombra como un rectángulo perfecto, que sería lo óptimo intuitivamente). A partir de cierto valor de $\sigma$ ($\sim 1$) obtenemos una forma de $S$, lo que captura muy bien la informaciónde de la superficie enfocada desde un lateral, pero no la extiende bien, que sería nuestro objetivo.
Por otro lado, utilizando KNN para computar la matriz de afinidad obtenemos un resultado relativamente parecido mientras aumentamos el número de vecinos utilizados. Sin embargo, las formas de $V$ obtenidas parecen extender peor la superficie a lo ancho para valores bajos de $k$. Además, si bien nos acercamos a obtener esa forma de $S$ como ocurre en RBF, esta es mucho menos consistente y la forma queda menos definida.
Hasta ahora hemos visto cómo trabajar con una matriz de datos que transformamos en un grafo para aplicar un método de reducción de dimensión espectral, pero todos los algoritmos que hemos vistos podrían aplicarse directamente a un grafo con cierta información. Para este ejercicio vamos a trabajar con el grafo del club de karate de Zacarías (un ejemplo clásico en teoría de grafos).
Para leer los datos necesitaréis instalar la librería igraph y cairo, para lo que podéis usar los comandos:
pip install python-igraph
pip install pycairo
Y si os da error porque faltan librerías deberéis instalar:
sudo apt install libcairo2-dev
sudo apt install python3-dev
Veamos primero que pinta tienen los datos:
import igraph
g = igraph.Graph.Read_GML("input/karate.gml")
igraph.summary(g)
# Información de una arista
print(g.es[0])
# Acceso a la información del primer nodo del grafo
print(g.vs[g.es[0].tuple[0]])
# Acceso a la información del segundo nodo del grafo
print(g.vs[g.es[0].tuple[1]])
IGRAPH U--- 34 78 --
+ attr: id (v)
igraph.Edge(<igraph.Graph object at 0x0000026836942E58>, 0, {})
igraph.Vertex(<igraph.Graph object at 0x0000026836942E58>, 0, {'id': 1.0})
igraph.Vertex(<igraph.Graph object at 0x0000026836942E58>, 1, {'id': 2.0})
El siguiente código permite visualizar el grafo dado. En este caso permitimos que se seleccione de manera automática el mejor algoritmo de visualización, etiquetamos los nodos por su identificador y lo coloreamos (todo del mismo color pues a priori no conocemos las comunidades).
layout = g.layout("auto")
g.vs["label"] = [int(idx) for idx in g.vs["id"]]
igraph.plot(g, layout=layout, bbox=(0, 0, 300, 300))
Puesto que los lados de este grafo no tienen pesos asociados (la forma más sencilla de apreciar esto es observar los datos en el documento proporcionado, karate.gml), la matriz de pesos del grafo será la matriz de adyacencia del mismo. Para obtenerla, recorremos los lados del grafo y y añadimos las aristas correspondientes a nuestra matriz.
W = np.array( g.get_adjacency().data )
print(W)
[[0 1 1 ... 1 0 0] [1 0 1 ... 0 0 0] [1 1 0 ... 0 1 0] ... [1 0 0 ... 0 1 1] [0 0 1 ... 1 0 1] [0 0 0 ... 1 1 0]]
Creamos una función que divide el grafo elegido (a partir de la matriz de pesos pasada por parámetro, W), en tantos clusters como elijamos (n_clusters) y proyectando a tantas dimensiones como queramos (n_components). Variando el parámetro method podemos elegir si utilizamos la implementación de Spectral Clustering de Sklearn o la nuestra.
Mostramos en primer lugar la nuestra y en segundo lugar, la de Sklearn.
def cluster_graph(W, n_clusters=2, n_components=2, method='ours'):
if method == 'ours':
sc = SpectralClustering(
n_clusters = n_clusters,
n_components = n_components,
affinity = 'precomputed',
precomputed_W = W,
)
y_pred = sc.fit_transform()
elif method == 'sklearn':
sc_sklearn = SpectralClusteringSklearn(
affinity = 'precomputed',
n_clusters = n_clusters,
n_components = n_components
)
y_pred = sc_sklearn.fit_predict(W)
else:
raise ValueError('Unknown method.')
# Swap labels 0 and 1 to deterministicly assign the same
# label to the first node
if y_pred[0] == 1:
y_pred[y_pred < 2] = (y_pred + 1) % 2
return y_pred
y_pred = cluster_graph(W, n_components=2)
g.vs['color'] = my_colors[y_pred]
igraph.plot(g, layout=layout, bbox=(0, 0, 300, 300))
y_pred = cluster_graph(W, method='sklearn')
g.vs['color'] = my_colors[y_pred]
igraph.plot(g, layout=layout, bbox=(0, 0, 300, 300))
Vemos como nuestra aproximación (reduciendo a dos dimensiones y tomando dos clusters) no captura tan bien como Sklearn la comunidad inferior derecha. Explorando con distinto número de dimensiones he encontrado que la mejor dimensión a la que proyectar para obtener ambas comunidades es n_components=1, siendo casi equivalente a a la de Sklearn pues sólo se diferencian en un nodo, el $20$.
Aunque puede parecer contraintuitivo proyectar en una única dimensión antes de realizar el clustering, esto es esencialmente proyectar en una recta y separar los dos clusters cortando la recta en algún punto.
y_pred = cluster_graph(W, n_components=1)
g.vs['color'] = my_colors[y_pred]
igraph.plot(g, layout=layout, bbox=(0, 0, 300, 300))